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INTRODUCTION 


During the third quarter progress continued along two fronts, 
the fluid dynamic model and the flame model of droplet combustion. 
This quarterly report describes some of the effort that is being 
generated in these two model building areas . 

The first section describes some numerical experiments 
carried out with the fluid dynamic pancake model of the combustor 
The experiments include both non-reacting and reacting flow con- 
ditions. The status of program COMB is also discussed. 

The second section describes recent work in producing a 
simplified droplet evaporation and combustion analysis. Some 
recent results towards completing the Peskin-Wise model relating 
to the modified flame surface analysis are presented. 

The authors would like to again acknowledge the creative 
programming effort that Harold Schechter has provided during the- 
course of this research program. 
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Fluid Dynamic Model 


A. Finite Amplitude Transverse Waves in a Cylindrical Chamber 

We have continued exploring and analyzing the nonlinear motion 
of finite amplitude waves in the transverse plane of a cylindrical 
channel. In the last quarterly report (number 2) we reported on 
a calculation which exhibited continuous behavior of the fluid 
properties for all time. The initial data is prescribed by rela- 
tions (traveling wave solution to the wave equation in cylindrical 
coordinates) which are presented on page 19 of quarterly report 
number 2. The initial pressure field and velocity field are shown 
in figures 1 and 2 and correspond to finite values for the ampli- 
tude parameter in the above equations. It should be noted that to 
obtain the dimensional value of pressure the nondimensional pressure 
is to be multiplied by yPco • this test case, however, the maxi- 

mum value of the pressure is 590 psia while the minimum value is 
10 psia (here p^ = 300 psia and y = 1.2). The numerical values 
shown in figures 1 and 2 correspond to this case. For the previous 
reported calculation the maximum (minimum) pressure was 450 (250) 
while p^ = 300 psia, and there the solution remained continuous 
for all time. 

As figure 3 shows, for the present calculation, a shock wave 
has formed after about 1.22 units of nondimensional time (time is 
normalized by R/a^ - .131 millisec) . The wave is strongest at the 
wall and decreases in strength as it curves inward into the chamber. 
Swirls seen behind the wave are due to numerical oscillations, a 
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feature inherent in these finite difference equations. Figure 4 
shows the strong discontinuous behavior of the velocity field in 
the neighborhood of the shock. Now figure 5, taken at 400 cycles 
of computation, is an intermediate stage of the flow field develop- 
ment. Two pressure peaks are evident at a nondimens ional radius 
of R = 1 of the combustion chamber. There is also a vast region 
of the .chamber where the pressure is approximately uniform and it 
is clear that figure 6 shows a strong induced velocity field which 
is predominently tangential in direction. The pressure field 
shown in figure 7 shows clearly a symmetric compression/rare faction 
wave with a peak pressure of about 1120 psia. The associated 
velocity field, shown in figure 8, is also. locally symmetric about 
the pressure wave. Velocities increase and then decrease uniformly 
in the direction of motion while the maximum velocity is at the 
boundary. 

One also notices that the velocity vectors on the leading edge 
of the wave point towards the wall while those at the trailing edge 
tend to point away from the wall. 

Figures 9 and 10 show essentially • the same picture as in the 
two previous plots; here the peak pressure is 1250 psia while the 
minimum pressure is 160 psia. The uniform portion of the flow 
field has a mean pressure level of about 220 psia. The wave has 
now undergone about one and two thirds rotations in 0.694 milli- 
seconds. Figures 11 through 14 again show that the wave has the 
same basic structure as seen before, however, the peak pressure is 
decreasing to 1018 psia. 
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Figures 15 and 16 show the pressure field with multiple 
waves, waves which then seem to coalesce and lead to the flow 
field given by the next three figures. 

It should be noted that the 'rev' counter value given at the 
bottom of the accompanying figures is computed from the relative 
position of the velocity vector at r = 0 from its starting posi- 
tion. For small values of the amplitude of the pressure wave, 
the vector rotates in phase with the wave. However, for large 
values of the amplitude, the fluid motion no longer is locked onto 
this convenient measure, but rather, there is in general a varia- 
tion in phase between the rotation of the wave and the velocity 
vector at r = 0 . At 2200 cycles (figure 18 corresponds to a 
physical time of 1.94 millisec) the wave has rotated about 4.8" 
times and the maximum amplitude has gone back up to 1350 psia. It 
is striking that the picture of the whole flow field again looks 
quite like that in figure 9. Apparently the fluid exhibits time 
wise dependent periodic behavior with waves of finite amplitude. 

The average period can be computed and is found to be 405 y-sec. 

Figures 20 and 21 give the dimensional pressure on the wall 
of the combustion chamber for selected values of time corresponding 
to some of the previous figures. The curves are labeled with three 
numbers. The first is the nondimensional, time; the second is the 
rotation counter enclosed in parenthesis; and the third is the 
computation cycle count. 
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B . Energy Release in a Fluid with Finite Amplitude Waves, 

Simple Energy Source 

We have carried out to completion the calculation started in 
the second quarterly report and described on page 26 of that re- 
port. 

The calculation consists of adding an energy source term to 
the energy equation. There is no associated mass source term in 
the continuity equation. The energy term is dependent on the local 
pressure via 

E = const j (p-1) if p<l 

v=l otherwise 

It is also assumed that a velocity effect can be included by 
allowing y to depend on the magnitude of the velocity |u| = (u^+v^)^ 
via 

y = 1 if juj z 6 
y = 0 otherwise 

We chose 6 to be large enough so that the energy release occurred 
within the region of the pressure wave and nowhere else. 

Finally, with the experience gained with the calculation des- 
cribed in section E of the second quarterly, we turn on the energy 
source slowly, i.e., the actual term used in the difference equation 
is o)fi where to is monitonically increasing and satisfies O^tosl. We 
take u to be the value of the 'rev' counter during the first rotation 
of the wave, it is taken as unity for all subsequent rotations. 

The calculation starts off exactly as in the last section. 
Figures 22 and 23 show the extent of the development of the wave 
after about 1.8 rotations of the wave (note the velocity vector 
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at r = 0, the rev counter, has just passed the starting position 
Q = -7r/2 for the second time) and it is clear that the structure 
is quite similar to the non-reactive case. However, the peak 
pressure is 4350 psia. The velocity field is also similar to the 
non-reactive case. The maximum velocity in the wave corresponds 
to 9800 ft/sec. Of course, as the calculation proceeds, the 
pressure amplitude continues to grow reaching a value approximately 
thirty times the reference pressure while the maximum velocity 
reaches a value six times the reference sound speed after just 
four hundred more cycles of calculation. 
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C . Energy Release in a Fluid with Finite Amplitude Waves 

(including the Axial Flux Calculation) 

It is clear from the last section that unless energy is allowed 
to leave the system convectively , the pressure (as well as the 
specific internal energy) will continue to increase raonotonically . 
Program COMB was then modified to include, in an approximate way, 
a value for the axial flux of energy . 

Several approximations were tried. The first method, and 
perhaps the most obvious, involved in essentially finding the total 
energy (and mass) contained in the r-0 plane for each cycle of 
computation. The excess due to accumulation is parceled equally 
so that at each mesh point there exists a contribution to outflow 
which will exceed the inflow by just the amount that is accumulated. 

It was found that this procedure seems to work for small amplitude 
waves where only slight spatial variations in the energy release 
exist. However, for finite amplitude waves which initiate locally 
large influxes of mass and energy (due to increased transport proc- 
esses on the droplet spray) this procedure fails to work. Some 
regions of .the flow field where little reaction occurs see, 
according to this procedure, a relatively large outflow of energy 
(or mass) . If this occurs for a sufficient number of time steps, 
the flow pressure (or density) can become negative locally and 
indeed this has been our observation. It is clear that a weighting of 
the outflow according to the inflow would be more desirable. 

This was done for the simple energy release model given in 
the last section; the results being given below. The energy that 
was input due to combustion was allowed to reside in the r, 6 
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plane for one time step and then it was immediately considered to 
contribute to the axial flux. The following table gives a quick 
comparison of the value of PmaxAPoo of three calculations. The 
first column is an inert calculation with no energy term, the 
second with the simple energy model including an energy flux in 
the axial direction, and finally the third calculation is with the 
simple energy model but no convective flux in the z direction, 
i.e., the energy is locked in. 

Wave Amplitude (Pmax/YP<„) 

Cycle Inert Case With Axial Flux Trapped Case 


0 

1.63 

1.63 

1.63 

600 

3.25 

3.73 

7.18 

1000 

3.21 

3.53 

19.57 

1200 

3.05 

3.45 

29.65 

1400 

2.84 

3.91 


The final row of data 

corresponds to 

approximately 3.27 

rotations of the 

wave and 

it is clear that 

by including the 


flux, the rapid growth of the wave is controlled while the energy 
source itself supplies energy to the wave. 

It was found that for values of ten times the constant in 
the energy law in section B, the solution could not be continued 
beyond 300 cycles. This value is consistent with a physically 
real energy release near the injector. It is clear that the eval- 
uation of the axial flux is quite critical and indeed, the problem 
really requires a third dimension so as to mathematically pose the 
calculation correctly. 
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D, 


Finite Amplitude Waves and Simplified Modified Godsave 
Analysis 


In this section we show the results of a tangential pressure wave 
interacting with a system of droplets (mean diameter = 75ji) which 
evaporates and burns according to the simplified Godsave analysis 
presented in Section II of this report. Hence, there is both an 
energy .source for .the energy equation as well as a mass source 
included in the continuity equatipn. However, there is no axial 
flux calculation included in this section. 

The calculation consisted in specifying the flow field at 
t = 0. It is similiar to that given in figures 1 and 2 except 
that the maximum (minimum) pressure is 550 psia (50 psia) . The 
calculation proceeds without energy addition for as long as it 
takes the wave to make one complete rotation. Then the calculation 
incorporates the relations developed in Section II of this report 
for the evaporation and combustion of the droplets . The energy 
and mass source is turned on completely just as the wave is 
starting its second period of rotation. Figures 24 and 25 give 
the fluid dynamic field just before the droplets start interacting 
(excuse the extraneous hash for this series of pictures which is 
due to a malfunctioning solenoid valve hooked onto the pen on the 
Calcomp Plotter) . For this series of pictures the pressure is 
normalized with p«> (rather than yP<» as before) . Figure 26 gives 
the dimensional pressure p(R,9) at t = 3.3259 units at the outer 
wall of the chamber (r = R) , again just before the energy is 
switched on. The wave is moving from right to left. 
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The next two figures, 27 and 28, show the initial transient 
as the energy is switched on. Due to large radial components, 
figure 28, the pressure pulse moves inward, but still undergoing 
rotation as the next four figures indicate (figures 29-32) . The 
next sequence of figures (figures 33-38) show the transition to 
the compression/rarefaction pulse traveling along the chamber 
wall (except that strong radial motion flattens the pulse against 
the wall) . This is seen from figure 39 to the final figure 44. 

As can be seen from figures 28 and 30, strong radial motion is 
generated by the large release of energy and persists as transverse 
motion is interacting with the tangential motion. The resultant 
'steady' or time periodic pattern would require additional cal- 
culation. 
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E . Program COMB Status 

As is shown by figure 26 , we have incorporated the ability 
to selectively plot the tangential pressure distribution at any 
value of radius and at any time. This feature adds flexibility to 
the program for ease in interpretation of the computed results. 

It is anticipated that the final report will contain some 
instantaneous streakline pictures of 'particles' moving in the 
r-6 plane. The resultant motion is shown plotted as separate 
pictures after, say, one complete rotation of the wave. It is 
hoped that correlations can be made with these computed results 
and photographic techniques to be attempted by the JPL rocket 
combustion group. 

The toroidal motor model is almost debugged. Tests are now 
being run to establish a steady state flow in the 6-z direction 
with given arbitrary initial data. A simple approach has been 
taken so as to simulate a converging diverging nozzle. That it 
is necessary to incorporate a supersonic boundary condition at 
the outflow plane of this toroidal geometry is evident. For a 
liquid propellant motor the only other required boundary condition 
is that the normal component of velocity shall vanish on the injec- 
tor face. Plotting routines have been written which will be used 
to see isobars and vector fields in this geometry. Some preliminary 
results have been obtained at this time but are too incomplete in 
nature to present. 
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II. Droplet Evaporation and Combustion Analyses* 


A. The Rate of the N2H4/N2O4 Reaction 


It will be recalled that in order to utilize the modified 
flame surface analysis of Peskin and Wise (Ref. 1 ) , some estimate 
is required of the overall (i.e,, 'global') reaction rate constant 
for the N2H4/N2O4 reaction (Ref. 2 ). That is, assuming the reaction 
proceeds according to 

N2H^ + NO 2 ^ 2^2 + 2 H 2 O (1) 


some estimate of k^ is required. (NO2 is considered the oxidant 
rather than N2O4 as a result of the high level of dissociation of 
the latter under most engine operating conditions.) 

This estimate can be obtained from the streamtube analysis 
developed in Ref. 2 , wherein the relevant equations and assumptions 
are detailed. For reference purposes, the streamtube equations 
are -summarized here as well: 

Energy Conservation 

™ = 0 , h = constant 

dt ' 


Species Conservation 

dYj ^ £i 
dt p 

Equation of State 

P = 

RT 


( 3 ) 


( 4 ) 


Figures in this section again start with the number 1 . 
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Auxiliary Equations 


= 5 Yihi 


M = 



hi = 






(5) 

( 6 ) 

(7) 


where h^ = enthalpy (sensible plus chemical) of species i 

= rate of appearance (or disappearance) of species i 
Yi = mass fraction of species i 
Mj^ = molecular weight of species i 
p = mass density of the gas mixture 
p = pressure of the mixture 
T == temperature of the mixture 

= reference enthalpy of species i 
Cp^ = average specific heat of species i 
Tr£p = reference temperature (1000°K) 
t = time 

Preliminary results were shown in Ref. 2. More complete 
results are presented in figure 1. Defining the ignition delay 
time (tj-p) qualitatively as the time required for the principal 
oxidation reactions to occur, figure 1 yields the anticipated 
result that tjj^j is a strong function of the initial mixture tern-- 
perature. On the other hand, these results indicate that in the 
range examined, tjj^ is nearly independent of the pressure. As 
can be seen, at the higher pressures, the second-order reactions 
(which become more significant at higher pressures) produce an 
initial partial oxidation which results in a slight temperature 
increase. However, a leveling-off period is seen to subsequently 
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occur, followed by the principal , first-order oxidation processes 
which drive the reaction to completion. 

As a result, the reaction may be said to be first-order, with 
ignition delay times (taken from figure 1) as shown in figure 2. 

In this latter figure, is the initial mixture temperature. The 
approximate straight line behavior of In tj^ when plotted against 
reciprocal initial temperature is in accordance with elementary 
chemical kinetic thoery which predicts a function of the form 

tiD = a -exp (E/RTj^) (8) 

where a is a constant and E is an overall activation energy. 

From figure 2, E was determined to be about 26,800 calories/mole. 

As an initial estimate, the pre-exponential factor in the 
overall reaction rate constant was chosen to be the average of 
the pre-exponential factors for the ^2^4 oxidation reactions as 
deduced by Sawyer (see Ref, 2, Table I) . From this and the result 
obtained above, the initial estimate of the overall reaction rate 
constant was 

kf = 1.2 X 10^° exp (-26,800/RT) (9} 

To check this estimate, results using eq. 9 were checked 
against results obtained using the complete streamtube analysis. 

The comparison was made at high pressure (20 atm) , where the 
first-order approximation and eq. 9 would be least accurate. The 
results are shown in figure 3, with the streamtube result labeled 
'exact'. The global model does not, of course, predict the. initial 
reactions or the leveling off period after the initial second-order 
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reactions have occurred, but the general behavior and the delay 
time appears to be reasonable. Although the result for this par- 
ticular case could have been improved by adjustment of the pre- 
exponential factor in eq. 9, this was felt to be unwarranted, and 
eq. 9 was selected as a suitably accurate overall reaction rate 
constant for inclusion in the modified flame surface analysis. 
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Under the assumptions inherent in the modified flame surface 
analysis, as detailed in Refs. 2 and 3, the species conservation 
and energy equations were shown to reduce to 



At \ = 

'5'' V aF y 


where a = dimensionless mass burning rate = mp/(4nrpk/c ) 

_o ^ 

mp = mass burning rate in stagnant surrounding 

r^ = fuel droplet radius 
k = thermal conductivity 
Ti = r/r^ 

ii = 

Vj^ = stoichiometric coefficient of species i 

Np,Nq = reaction orders with respect to fuel and oxidizer 

b = k.r^/D 
^ D 

D = binary diffusion coefficient 
t = (T-T^)/(Q/Cp) 

Q = heat of reaction 

= temperature in the ambient gas mixture 

On the basis of the conclusions drawn in the previous section, 
Nq = 0 and k^ is given by eg. 9. When eq. 10 is written first for 
the fuel and then for the oxidizer, and the resulting equations 
subtracted, an integrable ordinary differential equation results 
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subject to the following boundary conditions: 



When eq . 10 written for the fuel is added to eq . 11, another 
integrable ordinary differential equation results subject to these 
boundary ' conditions ; 


= L a 

dnj ri=l Q 
^ n-»-“ “ ® 

The resulting equations can be manipulated to yield the 
following expression for the dimensionless mass burning rate; 


a = In 


1 





(12) 


where L = latent heat of vaporization of the fuel 

tg = dimensionless droplet surface temperature 
Yp,s ~ mass fraction of fuel at the droplet surface 

We now introduce the additional assumption that the flame 
zone of finite thickness can be replaced by an inf initessimal one 
represented by a Driac-delta function located such that the fuel 
and oxidant react stoichiometrically . As a consequence, the right- 
hand-sides of eqs. 10 (written for the fuel) and 11 become 


64 



n-b Yp • 6 (n-n*) 

where the asterisks refer to the values at the hypothetical flame 
surface and 5 is the delta function. (Note that the fact that 
Np = 0 has been used in writing this term.) 

Peskin (Ref. 1) shows that a solution to equation 10 with 
this right-hand-side is given by generalized function theory as 


^ Y* n 


*0 forn>n* . . 

* ^ , a/n* ^ ^ 

1-e forn<ri* 


(13) 


Therefore, for n=n*, Yp=Yj, and eq. 13 becomes 


* 

Yp = 


(1-e-a/n*) 

1+ ~ n*^ (1-e-a/n*) 


(14) 


Now when eg. 13 is solved for n=lf Yp=Ypg, and eq. 14 is used, 
the result is 


b* *2 e-a (1-e^^^*) (1-e ^Z^^*) 

Ypg = ^ ^ +l-e-^ 

1+ n*2 (i-e-a/n*) 

a 


(15). 


From the difference between eq. 10 written for the fuel and 
the oxidizer, and using the appropriate boundary conditions, the 
following expression is obtained for n*: 


n 


* 




+ Yq («) 

e-a + io 



(16) 
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From the sum of equations 10 (written for the fuel) and 11, 
again subject to the appropriate boundary conditions, the following 
expression can be derived for t*: 


t* = (^FS ^s) (e~^/^*-l) - (l-e"^/n*) 


~a , 
e -1 


b^^ 72 


(1-e-a/^*) 

in which equation 14 was also used. 

Finally, using the result of the previous section, b* can 
be written 


(17) 


b* = fcL. {l,2 X 10^° • exp{-26,800/RT*)l (18) 

D 

Equations 12, 15, 16, 17, and 18 constitute five equations 
in the five unknowns a, Yps , n*, T*, and b* . (Unlike Peskin 
(Refs. 1, 4, 5), we are not interested here in species and tem- 
perature profiles surrounding the droplet, but in the mass burning 
rate a; as a consequence, the equations developed above are 
particularly suited to this purpose.) 

Results using this analysis have been obtained for the 
following values of the parameters; 

Yo(~) = 0.25 

D = 10“^ ft^ 
sec 

1q = 23 (see eq. 1) 

16 


Q was calculated from the enthalpies obtained from eq. 7; 
L was obtained from the Clausius-Clapeyron equation using the 
vapor pressure data in Ref. 6; in going from a to m°, k was 
obtained from Ref. 6 and Cp from Ref. 7. 
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Results are shown in figures 4 , 5 , .and 6 . In figure 4 , the 
results are compared with the previous results obtained from the 
Godsav analysis (Ref. 2) and the experimental results of Ref. 8. 

As can be seen, the current analysis yields results which are 
about 80% higher over most of the range of r^^ than the corresponding 
Godsav value. This is due, at least in part, to the values chosen 
for Yq(“) and D, as well as the differences inherent in the two 
approaches. It is interesting to note, however, that the line 
corresponding to T = 1035°K, p = 300 psia (Curve B) passes right 
through the experimental points of Ref. 8. Unfortunately, Ref. 8 
does not indicate the 'level of pressure at which these experiments 
were conducted, leaving this possible correlation unresolved. 
However, the linear dependence of the burning rate on droplet 
diameter is retained by the modified flame surface analysis, in 
accordance with the experimental results. 

Figure 5 contrasts the ambient temperature dependence of this 
analysis and the Godsav analysis for the particular conditions 
noted. It was pointed out in Ref. 2, and can be seen in figure 5, 
that the Godsav analysis yields a very weak, almost linear ambient 
temperature dependence. Most experimental work indicates a 
stronger, non-linear dependence. The modified flame surface 
analysis does yield this latter type of dependence. It may be 
justifiably said, however, that the differences between the two 
analyses ^ these particular conditions are not sufficiently great 
to warrant using the current analysis in the fluid dynamics pro- 
gram. In particular, no evidence of extinction can be seen down 
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to ambient temperatures near the boiling point of the fuel. If 
this same result is discovered at other values of the ambient 
pressure as well, this will constitute ample justification for 
the use of the simpler Godsav analysis . 

The effect of ambient temperature on the values of , n* 
and T* are shown in figure 6 . The value of Yps^l as T^ increases , 
in accordance with the assumption of the Godsav, diffusion-flame 
approximation . The hypothetical flame surface moves away 

from the droplet as Too increases; this stems from the fact that 
in order to sustain a stoichiometric combustion process, the flame 
location must move away from the flame surface where a larger 
quantity of fuel vapor (Yps) is being produced. Finally, T* ex- 
hibits a linear dependence on To,, 
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c. 


The Reduced Godsav Analysis 


Preliminary runs which coupled the complete Godsav analysis 
to the fluid dynamics program indicated that the running times 
were excessively long. (Roughly, five times slower than with the 
elementary energy term previously used.) As a consequence, the 
results obtained to date from the Godsav analysis (Ref. 2) were 
used to formulate the following simplified version of this analysis: 

The weak, almost linear, temperature dependence of the dimen- 
sionless mass burning rate can be approximately expressed as 

asoo = (3.73 x 10"5) (T„) + 1.855 (19) 

where ^300 is the dimensionless burning rate at a chamber pressure 
dependence is approximately given by 

_/P» \ 0.0108 

^00 ~(jUoJ ( 20 ) 

Once 3 p is computed, the burning rate in stationary surroundings 
is obtained from 

mp = 2 'irk dj^ ap ( 21 ) 

where dj^ is the drop diameter, and the thermal conductivity and 
specific heat are functions of the liquid drop temperature . 

The drop temperature, Tj^, which is the boiling temperature 
at the chamber pressure, can be obtained to a good approximation 
from the integrated Clausius-Clapeyron equation 
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In _2 = {8.7 X 10 ^) f 1 

2131 ITT76 


where p is in psia, Tj^ is in °R, and the critical point conditions 
are 2131 psia, 1176 °R. 

Knowing Tj|^, K is obtained from a linear curve fit of the 
data in Ref. 6; 

K = 1 r(-8.33 X 10“^) (Tt-460) + 0.2107 I (23) 

3600 L n j 

and Cp is obtained from Ref. 7; 

Cp = 0.138 + 0.527 X 10"^ Tj^ - 0.120 x 10“^ (24) 

Prom Ap, the burning rate in the convective environment is 
then obtained from (see Ref. 3, eq. 21); 

u 1/3 

■— = 1 + 0.276 (Re)^ (Pr)-^/-^ (25) 


Under the engine conditions of interest, the following con- 
stant values are reasonable: 

Pr^/^ =0.6 

Lp = 540 
^ X5m 


Lq = 178.2 


Q = 12,400 


where Lp is the latent heat of the fuel, is the latent heat of 
the oxidizer, and Q is the heat of reaction. The appropriate terms 
in the tj/-vector of the fluid dynamics program may now be obtained 
as shown in Ref. 3 (eq. 22) . Results obtained using this approxi- 
mate formulation are shown elsewhere in this report. 
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INTRODUCTION 


During the third quarter progress continued along two fronts, 
the fluid dynamic model and the flame model of droplet combustion. 
This quarterly report describes some of the effort that is being 
generated in these two model building areas. 

The first section describes some numerical experiments 
carried out with the fluid dynamic pancake model of the combustor. 
The experiments include both non-reacting and reacting flow con- 
ditions. The status of program COMB is 'also discussed. 

The second section describes recent work in producing a 
simplified droplet evaporation and combustion analysis. Some 
recent results towards completing the Peskin-Wise model relating 
to the modified flame surface analysis are presented. 

The authors would like to again acknowledge the creative 
programming effort that Harold Schechter has provided during the 
course of this research program. 
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I. 


Fluid Dynamic Model 


A. Finite Amplitude Transverse Waves in a Cylindrical Chamber 

We have continued exploring and analyzing the nonlinear motion 
of finite amplitude waves in the transverse plane of a cylindrical 
channel. In the last quarterly report (number 2) we reported on 
a calculation which exhibited continuous behavior of the fluid 
properties for all time. The initial data is prescribed by rela- 
tions (traveling wave solution to the wave equation in cylindrical 
coordinates) which are presented on page 19 of quarterly report 
number 2. The initial pressure field and velocity field are shown 
in figures 1 and 2 and correspond to finite values for the ampli- 
tude parameter in the above equations. It should be noted that to 
obtain the dimensional value of pressure the nondimensional pressure 
is to be multiplied by yPoo* 1^ this test case, however, the maxi- 
mum value of the pressure is 590 psia while the minimum value is 
10 psia (here p^ = 300 psia and y = 1.2). The numerical values 
shown in figures 1 and 2 correspond to this case. For the previous 
reported calculation the maximum (minimum) pressure was 450 (250) 
while p^ = 300 psia, and there the solution remained continuous 
for all time. 

As figure 3 shows, for the present calculation, a shock wave 
has formed after about 1.22 units of nondimensional time (time is 
normalized by R/a^ ~ .131 millisec) . The wave is strongest at the 
wall and decreases in strength as it curves inward into the chamber. 
Swirls seen behind the wave are due to numerical oscillations, a 



} 

feature inherent in these finite difference equations. Figure 4 
shows the strong discontinuous behavior of the velocity field in 
the neighborhood of the shock. Now figure 5, taken at 400 cycles 
of computation, is an intermediate stage of the flow field develop- 
ment. Two pressure peaks are evident at a‘ nondimens ional radius 
of R = 1 of the combustion chamber. There is also a vast region 
of the .chamber where the pressure is approximately uniform and it 
is clear that figure 6 shows a strong induced velocity field which 
is predominently tangential in direction. The pressure field 
shown in figure 7 shows clearly a symmetric compression/rarefaction 
wave with a peak pressure of about 1120 psia. The associated 
velocity field, shown in figure 8, is also locally symmetric about 
the pressure wave. Velocities increase and then decrease uniformly 
in the direction of motion while the maximum velocity is at the 
boundary. 

One also notices that the velocity vectors on the leading edge 
of the wave point towards the wall while those at the trailing edge 
tend to point away from the wall. 

Figures 9 and 10 show essentially - the same picture as in the 
two previous plots; here the peak pressure is 1250 psia while the 
minimum pressure is 160 psia. The uniform portion of the flow 
field has a mean pressure level of about 220 psia. The wave has 
now undergone about one and two thirds rotations in 0.694 milli- 
seconds. Figures 11 through 14 again show that the wave has the 
same basic structure as seen before, however, the peak pressure is 
decreasing to 1018 psia. 
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Figures 15 and 16 show the pressure field with multiple 
waves, waves which then seem to coalesce and lead to the flow 
field given by the next three figures. 

It should be noted that the ’rev' counter value given at the 
bottom of the accompanying figures is computed from the relative 
position of the velocity vector at r = 0 from its starting posi- 
tion. For small values of the amplitude of the pressure wave, 
the vector rotates in phase with the wave. However, for large 
values of the amplitude, the fluid motion no longer is locked onto 
this convenient measure, but rather, there is in general a varia- 
tion in phase between the rotation of the wave and the velocity 
vector at r = 0 . At 2200 cycles {figure 18 corresponds to a 
physical time of 1.94 millisec) the wave has rotated about 4.8 
times and the maximum amplitude has gone back up to 1350 psia. If 
is striking that the picture of the whole flow field again looks 
quite like that in figure 9 . Apparently the fluid exhibits time 
wise dependent periodic behavior with waves of finite amplitude. 

The average period can be computed and is found to be 405 u-sec. 

Figures 20 and 21 give the dimensional pressure on the wall 
of the combustion chamber for selected values of time corresponding 
to some of the previous figures. The curves are labeled with three 
numbers. The first is the nondimensional time; the second is the 
rotation counter enclosed in parenthesis; and the third is the 
computation cycle count. 
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B . Energy Release in a Fluid with Finite Amplitude Waves, 

Simple Energy Source 

We have carried out to completion the calculation started in 
the second quarterly report and described on page 26 of that- re- 
port. 

The calculation consists of adding an energy source term to 
the energy equation. There is no associated mass source term in 
the continuity equation. The energy term is dependent on the local 
pressure via 

E = const I (p-1) if p<l 

v=l otherwise 

It is also assumed that a velocity effect can be included by 
allowing y to depend on the magnitude of the velocity |u| = (u^+v^)^ 
via 

y = 1 if |u| 2 ; 6 
y = 0 otherwise 

We chose 6 to be large enough so that the energy release occurred 
within the region of the pressure wave and nowhere else. 

Finally, with the experience gained with the calculation des- 
cribed in section E of the second quarterly, we turn on the energy 
source slowly, i.e., the' actual term used in the difference equation 
is (jjfi where u is monitonically increasing and satisfies O^uil. We 
take to to be the value of the 'rev* counter during the first rotation 
of the wave, it is taken as unity for all subsequent rotations. 

The calculation starts off exactly as in the last section. 
Figures 22 and 23 show the extent of the development of the wave 
after about 1.8 rotations of the wave (note the velocity vector 
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at r = 0, the rev counter, has just passed the starting position 
e = — ir/2 for the second time) and it is clear that the structure 
is quite similar to the non-reactive case. However, the peak 
pressure is 4350 psia. The velocity field is also similar to the 
non-reactive case. The maximum velocity in the wave corresponds 
to 9800 ft/sec. Of course, as the calculation proceeds, the 
pressure amplitude continues to grow reaching a value approximately 
thirty times the reference pressure while the maximum velocity 
reaches a value six times the reference sound speed after just 
four hundred more cycles of calculation. 
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- • Energy Release in a Fluid with Finite Amplitude Waves 

(including the Axial Flux Calculation)" 

It is clear from the last section that unless energy is allowed 
to leave the system convectively , the pressure (as well as the 
specific internal energy) will continue to increase raonotonically . 
Program COMB was then modified to include, in an approximate way, 
a value for the axial flux of energy. 

Several approximations were tried. The first method, and 
perhaps the most obvious, involved in essentially finding the total 
energy (and mass) contained in the r-0 plane for each cycle of 
computation. The excess due to accumulation is parceled equally 
so that at each mesh point there exists a contribution to outflow 
which will exceed the inflow by just the amount that is accumulated. 

It was found that this procedure seems to work for small amplitude 
waves where only slight spatial variations in the. energy release 
exist. However, for finite amplitude waves which initiate locally 
large influxes of mass and energy (due to increased transport proc- 
esses on the droplet spray) this procedure fails to work. Some 
regions of the flow field where little reaction occurs see, 
according to this procedure, a relatively large outflow of energy 
(or mass). If this occurs for. a sufficient number of time steps, 
the flow pressure (or density) can become negative locally and 
indeed this has been our observation. It is clear that a weighting of 
the outflow according to the inflow would be more desirable. 

This was done for the simple energy release model given in 
the last section; the results being given below. The energy that 
was input due to combustion was allowed to reside in the r, 9 
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plane for one time step and then it was immediately considered to 
contribute to the axial flux. The following table gives a quick 
comparison of the value of Pjnax/YP« three calculations. The 
first column is an inert calculation with no energy term, the 
second with the simple energy model including' an energy flux in 
the axial direction, and finally the third calculation is with the 
simple energy model but no convective flux in the z direction, 
i.e., the energy is locked in. 

Wave Amplitude (PmaxAP^.) 

Cycle Inert Case With Axial Flux Trapped Case 


0 

1.63 

1.63 

1, 

.63 

600 

3.25 

3.73 

7, 

.18 

1000 

3.21 

3.53 

19, 

,57 

1200 

3.05 

3.45 

29. 

,65 

1400 

2.84 

3.91 




The final row of data corresponds to approximately 3.27 
rotations of the wave and it is clear that by including the axial 
flux, the rapid growth of the wave is controlled while the energy 
source itself supplies energy to the wave. 

It was found that for values of ten times the constant in 
the energy law in section B, the solution could not be continued 
beyond 300 cycles. This value is consistent with a physically 
real energy release near the injector. It is clear that the eval- ■ 
nation of the axial flux is quite critical and indeed, the problem 
really requires a third dimension so as to mathematically pose the 
calculation correctly. 
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D . Finite Amplitude Waves and Simplified Modified Godsave 

Analysis 

In this section we show the results of a tangential pressure wave 
interacting with a system of droplets (mean diameter = 75p) which 
evaporates and burns according to the simplified Godsave analysis 
presented in Section II of this report. Hence, there is both an 
energy source for the energy equation as well as a mass source 
included in the continuity equation. However, there is no axial, 
flux calculation included in this section. 

. The calculation consisted in specifying the flow field at 
t = 0. It is similiar to that given in figures 1 and 2 except 
that the maximum (minimum) pressure is 550 psia (50 psia) . The 
calculation proceeds without energy addition for as long as it 
takes the wave to make one complete rotation. Then the calculation 
incorporates the relations developed in Section II of this report 
for the evaporation and combustion of the droplets . The energy 
and mass source is turned on completely just as the wave is 
starting its second period of rotation. Figures 24 and 25 give 
the fluid dynamic field just before the droplets start interacting 
(excuse the extraneous hash for this series of pictures which is 
due to a malfunctioning solenoid valve hooked onto the pen on the 
Calcomp Plotter) . For this series of pictures the pressure is 
normalized with p» {rather than yp„ as before) . Figure 26 gives 
the dimensional pressure p(R,9) at t = 3.3259 units at the outer 
wall of the chamber (r = R) , again just before the energy is 
switched on. The wave is moving from right to left. 
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The next two figures, 27 and 28, show the initial transient 
as the energy is switched on. Due to large radial components, 
figure 28, the pressure pulse moves inward, but still undergoing 
rotation as the next four figures indicate {figures 29-32) . The 
next sequence of figures (figures 33-38) show the transition to 
the compression/rarefaction pulse traveling along the chamber 
wall (except that strong radial motion flattens the pulse against 
the wall). This is seen from figure 39 to the final figure 44. 

As can be seen from figures 28 and 30, strong radial motion is 
generated by the large release of energy and persists as transverse 
motion is interacting with the tangential motion. The resultant 
'steady' or time periodic pattern would require additional cal- ' 
culation. 
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E. 


Program COMB Status 


As is shown by figure 26, we have incorporated the ability 
to selectively plot the tangential pressure distribution at any 
value of radius and at any time. This feature adds flexibility to 
the program for ease in interpretation of the computed results. 

It is anticipated that the final report will contain some 
instantaneous streakline pictures of 'particles' moving in the 
r-6 plane. The resultant motion is shown plotted as separate 
pictures after, say, one complete rotation of the wave. It- is 
hoped that correlations can be made with these computed results 
and photographic techniques to be attempted by the JPL rocket 
combustion group. 

The toroidal motor model is almost debugged. Tests are now 
being run to establish a steady state flow in the 6-z direction 
with given arbitrary initial data. A simple approach has been 
taken so as to simulate a converging diverging nozzle. That it 
is necessary to incorporate a supersonic boundary condition at 
the outflow plane of this toroidal geometry is evident. For a 
liquid propellant motor the only other required boundary condition 
is that the normal component of velocity shall vanish on the injec- 
tor face. Plotting routines have been written which will be used 
to see isobars and vector fields in this geometry. Some preliminary 
results have been obtained at this time but are too incomplete in 
nature to present. 
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II. Droplet Evaporation and Combustion Analyses* 


A. The Rate of the N;>H4/N2Q4 Reaction 


It will be recalled that in order to utilize the modified 
flame surface analysis of Peskin and Wise (Ref. 1 ) , some estimate 
is required of the overall (i.e., 'global') reaction rate constant 
for the N2H4/N2O4 reaction (Ref. 2 ). That is, assuming the reaction 
proceeds according to 


N2H4 


+ NO2 



3 

2^2 + 2 H 2 O 


( 1 ) 


some estimate of kj is required. (NO2 is considered the oxidant 
rather than N2O4 as a result of the high level of dissociation of 
the latter under most engine operating conditions.) 

This estimate can be obtained from the streamtube analysis 
developed in Ref. 2 , wherein the relevant equations and assumptions 
are detailed. For reference purposes, the streamtube equations 
are summarized here as well: 

Energy Conservation 

^ = 0 , h = constant 

dt 


Species Conservation 

dYi _ ^i 
dt p 

Equation of State 

p = PM 
RT 


( 3 ) 


( 4 ) 


Figures in this section again start with the number 1 . 
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Auxiliary Equat ions 



(5) 

(6) 
( 7 ) 


where = enthalpy (sensible plus chemical) of species i 

r^ = rate of appearance (or disappearance) of species i 
= mass fraction of species i 
- molecular weight of species i 
p = mass density of the gas mixture 
p = pressure of the mixture 
T = temperature of the mixture 

= reference enthalpy of species i 
Cpi = average specific heat of species i 
"^REF ~ reference temperature (lOOO'^K) 
t = time 

Preliminary results were shown in Ref. 2. More complete 
results are presented in figure 1. Defining the ignition delay 
time (tjjj) qualitatively as the time required for the principal 
oxidation reactions to occur, figure 1 yields the anticipated 
result that tjj^ is a strong function of the initial mixture tem- 
perature. On the other hand, these results indicate that in the 
range examined, tjp is nearly independent of the pressure. As 
can be seen, at the higher pressures, the second-order reactions 
(which become more significant at higher pressures) produce an 
initial partial oxidation which results in a slight temperature 
increase. However, a leveling-off period is seen to subsequently 







occur., followed by the principal , first-order oxidation processes 
which drive the reaction to completion. 

As a result, the reaction may be said to be first-order, with 
ignition delay times (taken from figure 1) as shown in figure 2. 

In this latter figure, is the initial mixture temperature. The 
approximate straight line behavior of In tj^ when plotted against 
reciprocal initial temperature is in accordance with elementary 
chemical kinetic thoery which predicts a f' iction of the form . 

tiD = a.- exp (E/RTj^) (8) 

where a is a constant and E is an overall activation energy. 

From figure 2, E was determined to be about 26,800 calories/mole, 

As an initial estimate, the pre-exponential factor in the 
overall reaction rate constant was chosen to be the average of 
the pre-exponential factors for the N 2 H^ oxidation reactions as 
deduced by Sawyer (see Ref. 2, Table I). From this and the result 
obtained above, the initial estimate of the overall reaction rate 
constant was 

kf = 1.2 X 10^° exp (-26,800/RT) (9) 

To check this estimate, results using eq. 9 were checked 
against results obtained using' the complete streamtube analysis. 

The comparison was made, at high pressure (20 atm) , where the 
first-order approximation and eq. 9 would be least accurate. The 
results are shown in figure 3, .with the -streamtube result labeled 
'exact'. The global model does not, of course, predict the initial 
reactions -or the leveling off period after the initial second-order. 


59 















reactions have occurred, but- the general behavior and the delay 
time appears to be reasonable. Although the result for this par- 
ticular case could have been improved by adjustment of the pre- 
exponential factor in eg. 9, this was felt to be unwarranted, and 
eq. 9 was selected as a suitably accurate overall reaction rate 
constant for inclusion in the modified flame surface analysis. 
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Under the assumptions inherent in the ed flame surface 

analysis, as detailed in Refs. 2 and 3, the species conservation 
and energy equations were shown to reduce to 




.( 10 ) 


/ \ o ^o 

/n2 at \ b Yj, Yo 

'3" 'Si k ) <11> 

,o 

ffhere a = dimensionless mass burning rate = mp/(4irrQk/c ) 

r up 

O 

mp = mass burning rate in stagnant surrounding 
rj^ = fuel droplet radius 
k = thermal conductivity 
n = r/r^ 

±i = v^M^/VpMp 

= stoichiometric coefficient of species 1 

Np,Nq = reaction orders with respect to fuel and oxidizer 

b = k.r^/D 
^ D 

D = binary diffusion co.efficient 
t = (T-T„)/{Q/Cp) 

Q = heat of reaction 

= temperature in the ambient gas mixture 
On the basis of the conclusions drawn in the previous section, 
Nq = 0 and k£ is given by eq. 9. When eq . 10 is written first for 
the fuel and then for the oxidizer, and the resulting equations 
subtracted, an integrable ordinary differential equation results 
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subject to the following boundary conditions; 


§_ n- = a yJ 

du r^=l ^n=l 



When eq. 10 written for the fuel is added to eg. 11, another 
integrable ordinary differential equation results subject to these 
boundary conditions ; 


= L a 

dnj n=l Q 


Tl -»-00 — 0 

The resulting equations can be manipulated to yield the 
following expression for the dimensionless mass burning rate; 


a 




n 



(12) 


where L = latent heat of vaporization of the fuel 

tg = dimensionless droplet surface temperature 
Yp^s “ mass fraction of fuel at the droplet surface 

We now introduce the additional assumption that the flame 
zone of finite thickness can be replaced by an inf initessimal on 
represented by a Driac-delta function located such that the fuel 
and oxidant react stoichiometrically . As a consequence, the rig 
hand-sides of eqs. 10 (written for the fuel) and 11 become 
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^ b 


(S (n-ri*) 


where the asterisks refer to the values at the hypothetical flame 
surface, and 6 is the delta function; (Note that the fact that 
Nq = 0 has been used in writing this term.) 

Peskin- (Ref. 1) shows that a solution to equation 10 with 
this right-hand-side is given by generalized function theory as 

l-e^/*^ forn>n* 




Y* n*2 a 


1 Vn* ^ * 

1-e forn<n* 




(13)- 


Therefore, for n=n*, Yp=Xp and eq..l3 becomes 

yJ = (l-e~^/^*^ 

1+ — (l-e“V’i*) 

3 


(14) 


Now when eq. 13 is solved for ri=l, Yp=Ypg, and eq. 14 is used, 
the result is 


b* *2 

Y = a ^ 
^FS “ 


e-a d-e^/’^*) (l-e"^/"^*) 


+l-e 


-a 


1 + ^ n*2 (l-e-a/ri*) 

d. 


(15) 


From the difference between eq. 10 written for the fuel and 
the oxidizer, and using the appropriate boundary conditions, the 
following expression is obtained for n*: 


ip Ypg + Yq (oo) 

Yp (“) e-a + i^ Yp3 


(16). ■ 
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From the sum of equations 10 (written for the fuel) and 11, 
again subject to the appropriate boundary conditions, the following 
expression can be derived for t* : 


t* = (^FS ^s) (e~^/^*-l) - (i-e~^/n*j 


“ 3 - 1 
e -1 


(17) 


in which equation 14 was also used. 

Finally, using the result of the previous section, b* can 
be written 


b* = [[l.2 X 10^*^ • exp(-26,800/RT*)J (18)' 

D 

Equations 12, 15, 16, 17, and 18 constitute five equations 
in the five unknowns a, Ypg , n*, T* , and b* . (Unlike Peskin 
(Refs. 1, 4, 5) , we are not interested here in species and tem- 
perature profiles surrounding the droplet, but in the mass burning 
rate a; as a consequence, the equations developed above are 
particularly suited to this purpose.) 

Results using this analysis have been obtained for the 
following values of the parameters: 

Yo(«) = 0.25 

D = 10"^ ft^ 
sec 

i^ = 23 (see eq. 1) 

16 

Q was calculated from the enthalpies obtained from eq. 7j 

L was obtained from the Clausius-Clapeyron equation using the 

.0 , 

vapor pressure data in Ref. 6; in going from a to m^, k was 
obtained from Ref. 6 and Cp from Ref. 7. 
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Results are shown in figures 4, 5, and 6. In figure 4, the 
results are compared with the previous results obtained from the 
Godsav analysis (Ref. 2 ) and the experimental results of Ref, 8. 

As can be seen, the current analysis yields results which are 
about 80% higher over most of the range of r^^ than the corresponding 
Godsav value. This is due, at least in part, to the values chosen 
for Yo(~) and D, as well as the differences inherent in the two 
approaches. It is interesting to note, however, that the line 
corresponding to T = 1035*^K, p = 300 psia (Curve B) passes right 
through the experimental points of Ref. 8. Unfortunately Ref . 8 
does not indicate the level of pressure at which these experiments 
were conducted, leaving this possible correlation unresolved. 
However, the linear dependence of the burning rate on droplet 
diameter is retained by the modified flame surface analysis, in 
accordance with the experimental results. 

Figure 5 contrasts the ambient temperature dependence of this 
analysis and the Godsav analysis for the particular conditions 
noted. It was pointed out in Ref. 2, and can be seen in figure 5, 
that the Godsav analysis yields a very weak, almost linear ambient 
temperature dependence. Most experimental work indicates a ■ 
stronger, non-linear dependence. The modified flame surface 
analysis does yield this latter type of dependence. It may be 
justifiably said, however, that the differences between the two 
analyses at these particular conditions are not sufficiently great 
to warrant using the current analysis in the fluid dynamics pro- 
gram. In particular, no evidence of extinction can be seen down 
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to ambient temperatures near .the boiling point of the fuel. If 
this same result is discovered at other values of the. ambient 
pressure as well, this will constitute ample justification for 
the use of the simpler Godsav analysis. 

The effect of ambient temperature on the values of 
and T* are shown in figure 6. The value of increases, 

in accordance with the assumption of the Godsav, diffusion-flame 
approximation (^ps=l) * The hypothetical flame surface moves away- 
from the droplet as T„ increases; this' stems from the fact that 
in order to sustain a stoichiometric combustion process, the flame 
location must move away from the flame surface where a larger 
quantity of fuel vapor (^^3) is being produced. Finally, T* ex- 
hibits a linear dependence 
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c. 


The Reduced Godsav Analysis 


Preliminary runs which coupled the complete Godsav analysis 
to the fluid dynamics program indicated that the running times 
were excessively long. (Roughly, five times slower than with the 
elementary energy term previously used.) As a consequence, the 
results obtained to date from the Godsav analysis (Ref. 2) were 
used to formulate the following simplified version of this analysis; 

The weak, almost linear, temperature dependence of the dimen- 
sionless mass burning rate can be approximately expressed as 


a300 = (3-73 X lO'S) (T^) + 1.855 (19) 

where s^qq is the dimensionless burning rate at a chamber pressure 
dependence is approximately given by 


^300 



0.0108 


( 20 ) 


Once ap is computed, the burning rate in stationary surrounding; 
is obtained from 


riip = 2 irk Sp 

where dr is the drop diameter, and the thermal- conductivity and 
specific heat are functions- of the liquid drop temperature. 

The drop temperature, T^, which is the boiling temperature 
at the chamber pressure, can be obtained to a good approximation 
from the integrated Clausius-Clapeyron equation 
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In J 2 = (8.7 X 10-^) f 1 

2131 lll76 


where p is in psia, Tj^ is in °R, and the critical point conditions 
are 2131 psia, 1176 °R. 

Knowing Tj^, K is obtained from a linear curve fit of the 


data in Ref . 6 : 

K = 1 

3600 


|^(-8.33 X 10 (Tl-460) + 0.2107^ 


and c„ is obtained from Ref. 7: 

Ir 

c = 0.138 + 0.527 X 10"^ T, - 0.120 x 10"^ (24) 

P ^ 

From m°, the burning rate in the convective environment is 
F 

then obtained from (see Ref. 3, eq. 21) : 

1 + 0.276 (Re)^ (Pr) (25) 

mp 

Under the engine conditions of interest, the following con- 
stant values are reasonable: 

= 0.6 

% = 540 ^ 


Lq = 178.2 


Q = 12,400 _ 

Ibm 

where Lp is the latent heat of the fuel, Lq is the latent heat of 
the oxidizer, and Q is the heat of reaction. The appropriate terms 
in the iJ;-vector of the fluid dynamics program may now be obtained 
as shown in Ref. 3 (eq. 22) . Results obtained using this approxi- 
mate formulation are shown elsewhere in this report. 
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